function [H, er] = SWsymnmf(M, c)
% Smallwhite Symmetric NMF Package (swsymnmf)
% by Sun Sibai (Smallwhite <niasw@pku.edu.cn>)
%
% This is the entrance function of Symmetric NMF. The algorithm is ANLS.
%
% Variables:
%   >> M: sparse symmetric matrix to be factorized (Similarity Matrix)
%   >> c: number of dimensions of factors (Cluster Number)
%   << H: factor matrix (Cluster Result)
%   << er: residual error (to be minimized)
%
% References:
%   * symnmf_anls [http://www.cc.gatech.edu/~dkuang3/]
%   > by `Da Kuang`, `Chris Ding`, `Haesun Park`,
%   > Symmetric Nonnegative Matrix Factorization for Graph Clustering,
%   > The 12th SIAM International Conference on Data Mining (SDM '12), pp. 106-117.
%   > (The algorithm is from here. But that one does not support sparse matrix. So I have to rewrite them.)
%
% Basic Mechanism:
%   min_{H,W} er(H) = ||M - WH'||_F^2 + a * ||W-H||_F^2
%   H_{j,k}>=0, a = max(M_{j,k})
%
% New Features:
%   * normalize H on each iteration: if no normalization is applied, H will decrease to 0 which may cause clustering fail.
%
 MAXITERLIM = 10000; % limit max iteration to prevent dead loops

 n = size(M,1);
 if (n~=size(M,2))
  error('(ToT) Input Error: M should be symmetric.');
 end
 if (c<=0 || c>n || c~=int64(c))
  error('(ToT) Input Error: c should be a positive integer no more than dim(M).');
 end

 disp(sprintf('(^u^) dim(M) = %d x %d',n,n));
 disp(sprintf('(^u^) dim(H) = %d x %d',n,c));
 [it1,it2,val] = find(M); % retrieve all nonzero elements of sparse matrix M
 a = max(val); % use max(M) as the parameter a

 if (a<=0)
  error('(ToT) Input Error: M should contain positive elements.');
 end

 H = 2*sqrt(mean(val)/c) * rand(n,c); % initialization as 2 sqrt(avg(M)/c) Random(n,c) [not sparse]
 W = H;
 Ic = a * speye(c); % c x c identity matrix

 % shorten updating time cost
 L = H'*H; R = M*H;

 for it = 1:MAXITERLIM

  W = nnlsm_blockpivot(L + Ic, (R + a * H)', 1, W')';
  L = W' * W;
  R = A * W;
  H = nnlsm_blockpivot(L + Ic, (R + a * W)', 1, H')';
  tempW = sum(W, 2);
  tempH = sum(H, 2);
  temp = a * (H-W);
  gradH = H * L - R + temp;
  L = H' * H;
  R = A * H;
  gradW = W * L - R - temp;

  er = swfnorm(M,H);

  disp(sprinf('(*o*) iter: %d  err: %g',it,er));

  if (it == 1)
   initgrad = sqrt(norm(gradW(gradW<=0|W>0))^2 + norm(gradH(gradH<=0|H>0))^2);
   continue;
  else
   projnorm = sqrt(norm(gradW(gradW<=0|W>0))^2 + norm(gradH(gradH<=0|H>0))^2);
   if projnorm < 1e-8 * initgrad
    break;
   end
  end

 end
end
